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The nonlinear development of the strong Buneman instability and the associated 
fast electron heating in thin current layers with Q e /u pe < 1 are explored. Phase 
mixing of the electrons in wave potential troughs and a rapid increase in temperature 
are observed during the saturation of the instability. We show that the motion of 
trapped electrons can be described using a Hamiltonian formalism in the adiabatic 
approximation. The process of separatrix crossing as electrons are trapped and de- 
trapped is irreversible and guarantees that the resulting electron energy gain is a true 
heating process. 
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The exploration of how waves and particles interact in strong turbulence has been an 
important challenge in plasma physics. 1-9 Using particle- in-cell simulations we explore the 
nonlinear development and nonlinear wave-particle interactions of the Buneman instability 
to reveal how particle acceleration and heating take place. The Buneman instability— is 
driven by the relative drift between ions and electrons. Its quasi-linear theory has been 
well understood, but strong Buneman turbulence is poorly understood though it has been 
widely discussed^— . The previous work largely focussed on linear kinetic diffusion theory 
with an approximation for strong turbulence which assumes the electron heating is caused 
by wave-particle resonance with weakly non-linear waves. Instead, we have found that 
rapid electron heating takes place near the saturation of the Buneman instability when 
the particle's bounce rate in the wave potential is far larger than the growth rate of the 
instability. As a consequence, the trapped particle's motion is approximately adiabatic. The 
heating is a consequence of coherent trapping, phase mixing and de-trapping of the particles. 
Our simulations also demonstrate the difference between the nonlinear development of the 
Buneman instability and an idealized adiabatically-growing wave. 

Electron heating as a result of the Buneman instability is associated with the intense 
electron current layers formed during magnetic reconnection^"— , shocks^ 1 ^ 1 ^ and turbu- 
lent energy cascades to sub-proton scales^ 9 - 1 ^. In particular, understanding how kinetic 
turbulence transfers momentum and energy is important for revealing the role of anomalous 
resistivity in magnetic reconnection, which has been a long-standing puzzle^ 1 ^. 

We propose a new mechanism that is responsible for extremely fast electron heating, in 
a few tens of electron plasma periods, during the nonlinear evolution of the Buneman insta- 
bility. The dynamics is dominated by the coherent trapping and de-trapping of streaming 
electrons (with drift Vd e ) in the nearly non-propagating electric field from the instability. The 
wave amplitude grows until nearly all of the streaming electrons have been trapped. Thus 
the electrostatic potential at saturation is approximately given by ecfi ~ m e v\ e . The bounce 
frequency Uj, = k ^Je(j)/m e ~ uj pe of electrons trapped in the potential <fi greatly exceeds 
the characteristic growth rate of the wave 7 ~ [m e /mi) l ^ijjp e . As a result, the electrons 
trapped in the growing potential behave adiabatically, preserving their phase space area 
as the wave amplitude slowly changes in time. Phase mixing of the electrons in the wave 
potential troughs guarantees that, as the wave amplitude decreases following saturation, the 
de-trapping of electrons leaves a distribution of particles that forms a velocity-space plateau 
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over the interval (— v^, Vde)- The process of separatrix crossing as electrons are trapped and 
de-trapped is irreversible and guarantees that the resulting electron energy gain ~ m e v\J1 
is a true heating process. We carry out 3D PIC simulations with strong electron drifts in 
an inhomogeneous current- carrying plasma with a guide field. We apply no external per- 
turbations to initiate reconnection, and consequently reconnection does not develop during 
the simulation. We specify the initial magnetic field to be B X /B Q = tanh[(y — L y /2)/wo], 
where -Bo is the asymptotic amplitude of B x , and wq and L y are the half-width of the initial 
current sheet and the box size in the y direction, respectively. The guide field B 2 = B 2 — B 2 
is chosen so that the total field B is constant. In our simulation, we take wq = O.ldi and 
B = \/26Bq, where di = c/u P i and u P i is the ion plasma frequency. The initial tempera- 
ture is T e0 = Tj = 0.04mjC^ , the mass ratio is mj/m e = 100, and fl e ~ 0.625w pe , where 
cao = Bq/ {Annomi) 1 / 2 is the asymptotic ion Alfven wave speed. The simulation domain 
has dimensions L x x L y x L z = 0.5 x 1 x 4 di with periodic boundaries in x and z and 
a conducting boundary in y. The initial electron drift along z, v&, ~ lOc^o, is above the 
threshold for triggering the Buneman instability. 

During the simulation the Buneman instability onsets around Q e t ~ 25.5 with a wave 
vector that is aligned along the magnetic field. The instability reaches its peak around 
Q e t ~ 40 and ceases around Q e t ~ 102, as indicated by the turbulence strength Var(E z ) 
in FigJTk (blue dot-dashed line). The electric field E z parallel to B z abruptly increases 
from a few E = caqBq/c to E z ~ 40 — 60E at Q e t ~ 40 and then falls to a value 
20E at Q e t ~ 64 (Fig{Tb). At the same time, the average parallel component of the 
electron temperature, {T ezz ) sharply increases, from 0.04 to 0.5, by more than a factor of 
10 while the ion temperature increases only slightly. () denotes an average over the mid- 
plane of the current sheet at y — L y /2. The electron drift velocity decreases from 9cao to 
~ Qcao- It is noteworthy in Fig. [TJi that the increase of (T ezz ) nearly matches the damping 
rate of the electron parallel kinetic energy W, which implies that the streaming kinetic 
energy is nearly fully converted into thermal energy, i.e. -^A(t>J e ) w -A(T ezz ), where the 
Boltzmann constant has been absorbed into T. Panel b in Fig. [TJ shows the electron velocity 
distribution function in the current sheet at Q e t = 25.5, 38.3, 51, 63.7 and 102. The narrow 
ion distribution function at Q e t = 102 is shown with a solid line. We can see that the 
electron velocity distribution functions f(v ez ) become flatter and broader at late times, but 
the significant change takes place during Q e t ~ 38 — 51. The electron distribution functions 
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FIG. 1. a: The time evolution of d(T ezz ) Idt (black solid line), \dW/dt\ = ^dijl JnT)/dt (black 

dot-dot-dashed line), and the parallel electric field variance Var(E z ) = (y 7 (E z — (E z )) 2 ) (blue dot- 
dashed line), where n e is the electron density and j ez is the z-component of the electron current 
density, b: Electron velocity distribution functions in the current sheet f{v ez ) at Q e t = 25.5, 38.2, 
51, 63.7 and 102 plotted as black dot-dashed, red short-dashed, orange long-dashed, blue dot-dot- 
dashed and black solid lines respectively. The narrow ion velocity distribution function f(vi z ) at 
£l e t = 102 is reduced by 5 times, c: Electric field E z parallel to B at £l e t =25.5 (blue dot-dashed), 
38.2 (black solid), 63.7 (red dashed). 

reveal that a few electrons are accelerated to very high velocity, which is a consequence 
of the inductive electric field E z that maintains the integrated current. The ion velocity 
distribution function is slightly affected by the turbulence. 

To fully understand the fast electron thermalization, we plot the electron distribution 
function in phase space (z, v ez ) at times Q e t = 25.5, 38.2, 51, and 63.7 in Fig. [2j We see that 
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FIG. 2. Electron distribution functions in phase space (z,v ez ) at Q. e t =25.5, 38.2, 51 and 63.7. 



at Q e t = 38.2, electron begin their trapped orbits in the potential wells and by Q e t = 51 they 
have phase mixed along their trapped orbits. The localized electric field structure remains 
intact. Phase mixing does not change significantly after Q e t = 63.7. The period of fast phase 
mixing is coincident with that of the flattening of electron distribution functions shown in 
Fig. [TJ The phase mixing occurs near the time of saturation when the change in the electric 
field is small. We now explore the physical mechanism that produces the fast phase mixing 
and electron heating. 

To reveal the physics behind this phase mixing and electron heating, we need to identify 
the properties of the instability driving the turbulence. We use a double drifting-Maxwellian 
kinetic model to trace the evolution of the instabilities in the simulation in the following 
way^ 1 ^ 1 ^. We fit the ion and electron distribution functions averaged over z at x = 0.025 in 
the mid-plane x — z of the current sheet at Q e t = 25.5, 38.2, 51, and 63.7, and then substitute 
the fits into the local dispersion relation derived from a double drifting-Maxwellian kinetic 
model for waves with fL <C u <C Q e : 
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[1 + J (A)e- A C e2 Z(Ce2)] = 0, (1) 



where Q = {u - k z v di )/kv ti , Qi = (w - k z v del )/k z v ztel , ( e2 = (u - k z v de2 )/k z v zte2 , X 



klvl te /2Ql, Z is the plasma dispersion function and 1$ is the modified Bessel function of the 
first kind with order zero. The thermal velocity of species j is defined by = 2T t j/rrij and 
the drift speed by which is parallel to the z direction. The electron temperature takes 
on different values parallel and across the magnetic field while the ions are assumed to be 
isotropic. 5 is the weight of the high velocity drifting Maxwellian. 

We numerically solve the dispersion relation in Eq. [1] and obtain the unstable modes at 
Q e t = 25.5, 38.2, 51, 63.7 and 102. We find that the Buneman instability dominates. The 
growth rate of the fastest growing Buneman mode decreases with time from 7 ~ 0.12u pe 
(close to the classic value 7 ~ (m c /2m i ) 1 / 3 u pc ~ 0.17u pe ) at Vt e t = 25.5 to 7 ~ 0.0Qu pe at 
Q e t = 38.2 and 7 ~ 0.006cj pe at Q e t = 102. The frequency of the fastest growing mode is 
about ujq ~ 0.013o;p e and its wavenumber k di decreases from 90 to 75. The phase speed 
v p increases slightly with time and has a value of v p ~ 0.05cao- A transient two stream 
instability with growth rate 7 ~ 0.006w pe appears at Q e t = 51 and is stable by Q e t = 63.7. 
An oblique lower hybrid instability develops with growth rate 7 ~ 0.02u pe after Q e t = 51. 

It's interesting to notice that during Q e t ~ 38 — 64, the typical parallel electric field 
is about 40^0 and the wavenumber of the fastest mode is k$di ~ 90. The corresponding 
bounce frequency is ujb = k Vb/V2 ~ oo pe , where ~ ^/2e<f)/m e ~ lOc^o- The bounce rate 
is more than ten times larger than the growth rate. Thus the amplitude of the electric field 
E z evolves slowly compared with oj^ 1 during this interval. By assuming the slowly evolving 
and slowly propagating (i.e. z — v p t ~ z) wave potential is <p(z, t), the electron Hamiltonian 
can be approximated as : 

H^vl-e^t). (2) 

Eq. (T5]) shows that during fl e t = 38 — 64, it is possible to choose z and v ez as two 
approximate Hamiltonian canonical coordinates so that the area S = j- § v ez dz enclosed 
by the electron's trajectory in phase space (z, v ez ) is an adiabatic invariant for trapped 
electrons, where v ez = ^j2(W{t, z) + e0)/m e and W(t, z) is the electron's total energy. 
With the slow variation of the electric field, the electron's trajectory in phase space (z,v ez ) 
becomes narrower in z and longer in v ez as E z increases and becomes wider in z and shorter 
in v ez as E z decreases. The electrons are trapped when the electric field grows and are de- 
trapped when the electric field decays. The trapping and de-trapping are non-adiabatic due 
to the change of the phase area inside and outside of the wave potential^. The final electron 
velocity depends on whether it crosses the upper or lower separatrix as it is de-trapped. The 
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FIG. 3. Panel (a) and (b) are electron distributions in phase space z,v ez at the middle and the 
end time of the test particle simulation with E max ~ 80. The black dashed lines in Panel (c) and 
(d) are the initial test electron velocity distribution functions. Panel (c) is at the middle time of 
the test particle simulations and panel (d) is at the end time of the simulations. The black solid 
lines are for the test particle simulation with E max ~ 80 and the red dot-dashed line are for the 
test particle simulation with E max ~ 40. 



upper (lower) separatrix crossing results in a positive (negative) velocity in the wave frame. 

To investigate how the adiabatic process converts kinetic energy into thermal energy 
through non-adiabatic separatrix crossings of the wave potential, we perform two test particle 
simulations with 5000 electrons in one single standing wave E = E^e 11 sin kz. We take 
kdi ~ 80 ~ ho] E z q and |7| are constant and small. They satisfy oob/\l\ ~ 80 at the peak 
value of E z . 7 > for the first half of the total duration and 7 < for the second half 
so that E grows and decays sufficiently slowly to guarantee that the motion of the trapped 
particles is adiabatic during the entire duration. The duration is determined by the peak 
value of E max . We investigate the cases with E max ~ 40 and E max ~ 80. The initial electron 
velocity distribution is a Maxwellian with a drift Vd e ~ 9cao and T e = 0.04mjC^ and the 
electron density is uniform in space. The value E max = 40 is similar to the peak value of E z 
observed when the PIC simulations can trap electrons with velocity below i^ e . E max = 80 
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is higher than the peak value of E z observed in the PIC simulations. The test single wave 
with E max = 80 can trap almost all of the electrons. The results are shown in Figj3j 

More and more electrons are trapped as the electric field slowly increases and the most 
are trapped at the peak value of E. The slight energy difference between two trapped 
electrons leads to a large separation in their phase angle around their trapped orbit since 
their angular velocity depends on energy. Thus, at the time of the maximum trapping 
the trapped electrons are nearly uniformly distributed along their trajectories as shown in 
panel a and the velocity distribution of trapped electrons become flat as shown in panel 
c. As E z decreases, the electron energy gain during trapping reverses and the electrons are 
eventually de-trapped with the same value of W. The total energy W is symmetric with 
respect to positive and negative velocity. Therefore, at the end of the simulation, due to 
the same probability for de-trapping at the positive and negative velocity (Figj3] panel b), 
a dip appears near zero velocity in the velocity distribution function shown in panel d. The 
red lines in panel c and d are for E max ~ 40 where E max is not strong enough to trap all of 
the electrons. As a result, the distribution functions are not completely symmetric around 
v ez = 0. The red line in panel c is similar to the electron velocity distribution function at 
the saturation stage of the Buneman instability displayed in FigJTb. 

However, even at Q e t = 102 in Fig. [Tb, the electron velocity distribution function keeps a 
similar shape to that at the saturation stage rather than show a dip near zero as seen in the 
test particle simulation. There are two reasons for the missing dip in the PIC simulations. 
First, the wave amplitude is not spatially uniform. As a result, the trapping/de-trapping is 
more complex than in the sample model. In Fig. HI we plot the constant energy contours of 
electrons in phase space at Q e t = 40 and 102 which correspond, respectively, to the peak of 
the Buneman instability and late time of turbulence. We can see that the islands between 
z G [0, 0.2] and z G [0.3, 0.4] are longer than those between z G [0.2, 0.3] and z G [0.4, 0.5] at 
Qit = 40 due to the corresponding variation of the electric field E z in z as shown in Fig. [lb. 
The long islands in phase space correspond to weak electric field and weak trapping; Second, 
after Q e t = 51, the electron two stream and Buneman instabilities remain unstable, albeit 
weaker, and trapping continues. At the late stage of the turbulence Q e t = 102, the islands 
in phase space (Fig. HJ) still exist even though the islands become longer. 

In the PIC simulation, the electron trapping velocity cover the range [—10, 10]cao ~ 
[—v de ,Vde\- The heating stops when the source of kinetic energy is completely drained, 
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FIG. 4. The contours of total electron energy in phase space (z,v ez ) at Q e t = 40, 102 from the 
PIC simulations. 

i.e. the electron distribution with velocity below become flat, as shown in Fig. [T] b. 
Nearly half of the kinetic energy is dissipated, and the final electron temperature can be 
estimated by AT ezz ~ AW ~ m e v^ e /2 ~ 0.4 which is consistent with what is observed in 
the simulation (Fig. [T^,). 

In reconnection, current sheets shrink as reconnection evolves and the Buneman insta- 
bility might occur in a wider current sheet with a reduced drift but a similar growth rate 
~ (m e /2m i ) 1 / 3 uj pe . The implications of these results in reconnection are being explored. 
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